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We  present  a  molecular  dynamics  computer  simulation  method  for  calculating 
equilibrium  constants  for  the  formation  of  physical  clusters  of  molecules. 

The  method  is  based  on  Hill's  formal  theory  of  physical  clusters.  In  the 
method,  a  molecular  dynamics  calculation  is  used  to  calculate  the  average 
potential  energy  of  a  cluster  of  molecules  as  a  function  of  temperature,  and 
the  equilibrium  constants  are  calculated  from  the  integral  of  the  energy  with 
respect  to  reciprocal  temperature.  The  method  is  illustrated  by  calculations' 
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I.  Introduction 

The  thermodynamic  stability  of  clusters  of  small 
numbers  of  atoms  or  molecules  can  be  described  quantita¬ 
tively  by  the  equilibrium  constants  for  the  formation  of 
the  clusters.  For  the  reaction 

n  A  -  A  , 
n* 

the  equilibrium  constant  is 

K  -  (A  )/ (A)n 
n  n 

where  parentheses  about  a  species  name  denotes  the  number 
density  of  the  species.  A  knowledge  of  these  equilibrium 
constants  is  important  in  a  number  of  areas  of  chemical 
physics.  First,  the  values  of  Kg  and  are  related  to  the 
second  and  third  virial  coefficients  of  the  substance. 

Second,  the  theory  of  nucleation  of  the  liquid  phase  (or 
the  solid  phase)  in  a  supersaturated  vapor  uses  these  equi¬ 
librium  constants  and  their  closely  related  free  energies. 
Third,  it  has  been  suggested  that  the  variable  continuous 
absorption  of  infrared  radiation  by  the  atmosphere  arises 
from  these  clusters.^  A  reliable  and  practical  method  for 
the  theoretical  calculation  of  cluster  formation  equilibrium 
constants  would  have  consequences  for  each  of  these  three  areas 
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The  calculation  of  the  equilibrium  population  distribution 

of  water  clusters  is  of  particular  current  importance  because 

2 

of  the  controversy  surrounding  the  recent  hypothesis  that 
water  clusters  are  much  more  prevalent  in  the  atmosphere 
than  had  been  assumed. 

Equilibrium  constants  for  reactions  such  as 

n  A  +  B  =  BA  , 
n 

where  A  is  water  and  B  is  an  ion  can  be  measured  by  mass 

3 

spectrometric  techniques.  A  theory  for  the  calculation  of 
such  constants  would  enable  the  experiments  to  be  interpreted 
quantitatively  in  terms  of  the  basic  interactions  between 
ions  and  water  and  between  water  molecules. 

The  classical  way  to  evaluate  equilibrium  constants  for 
clusters  of  identical  molecules  treats  the  cluster  as  a 
droplet  whose  free  energy  contains  bulk  and  surface  contri¬ 
butions.^’^  Such  a  theory  makes  a  connection  with  the  macro¬ 
scopic  thermodynamic  properties  of  the  material  but  does  not 
address  the  relationship  between  intermolecular  interactions 
and  the  stability  of  clusters. 

The  partition  functions  for  clusters  of  atoms  and 
molecules  can  in  principle  be  calculated,  leading  to  a  sta¬ 
tistical  mechanical  theory  of  cluster  formation  constants 
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that  is  similar  to  the  usual  statistical  mechanical  theory 
of  equilibrium  constants  for  chemical  reactions.  Hill's 
physical  cluster  theory^  provides  a  rigorous  formalism  for 
defining  the  cluster  equilibrium  constants  and  relating 
them  to  the  thermodynamics  of  a  gas. 

If  one  applies  to  clusters  the  same  approximations  that 
are  used  in  calculating  partition  functions  for  molecules 
in  the  gas  phase,  such  as  rigid  rotor  and  harmonic  oscillator 
approximations,  then  partition  functions  and  their  associated 
free  energies  and  equilibrium  constants  can  be  calculated. 
Such  calculations  are  most  appropriate  for  clusters  if  the 
temperature  is  low,  if  the  configuration  is  solid-like,  and 
if  there  is  a  unique  lowest  energy  configuration  that  domi¬ 
nates  the  thermodynamic  properties  at  low  temperatures. 

(See  Ref.  7  for  a  discussion  of  the  single  configuration 
approximation. ) 

Computer  simulation  methods,  such  as  the  molecular 
dynamics  and  Monte  Carlo  methods,  can  be  used  to  eliminate 
statistical  mechanical  approximations.  Many  studies  of 
atomic  and  molecular  clusters  have  been  made  using  these 
methods.  A  fundamental  problem  associated  with  using 

simulation  methods  for  the  evaluation  of  cluster  properties 
is  that  the  simulation  of  a  cluster  at  a  particular  tempera¬ 
ture  cannot  be  used  to  evaluate  the  equilibrium  constant 
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or  cluster  free  energy  at  that  temperature.  The  quantity 
that  is  most  easily  calculated  is  the  average  energy  of  a 
cluster  as  a  function  of  temperature.  Most  of  the  computer 
studies  listed  above  were  in  fact  calculations  of  the  energy 
and  not  the  free  energy.  The  energy  is  related  to  the 
temperature  derivative  of  the  free  energy,  and  so  calcula¬ 
tion  of  the  energy  for  a  range  of  temperatures  can  lead  to 
free  energies  by  numerical  integration,  but  the  results  con¬ 
tain  an  unknown  constant  of  integration  that  represents  the 
free  energy  in  the  state  at  which  the  integration  is  started. 

There  are  several  ways  around  this  difficulty.  One  is 

to  start  the  integration  at  a  very  low  temperature  state  and 

use  the  harmonic  approximation  and  the  single  configuration 

approximation  to  evaluate  the  free  energy  of  that  state.  A 

9 

second,  due  to  Lee  et  al . ,  makes  use  of  the  arbitrariness 
in  the  definition  of  a  cluster.  They  developed  a  way  of 
calculating  the  derivative  of  the  free  energy  with  respect 
to  the  radius  of  the  hypothetical  containment  sphere  used 
to  define  the  cluster.  For  infinite  containment  spheres, 
the  atoms  in  the  cluster  are  mostly  far  from  each  other  and 
analytical  methods  can  be  used  to  calculate  the  free  energy. 
Then  by  numerical  integration  with  respect  to  the  sphere 
radius,  they  are  able  to  find  the  absolute  free  energy  of  a 
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cluster  defined  with  a  finite  sphere  containment  radius. 

A  third  method  was  developed  by  Mruzik  e£  al.^  They  used 
a  coupling  parameter  method  and  evaluated  the  average 
potential  energy  of  interaction  of  one  molecule  with  the 
n  -  1  other  molecules  in  the  cluster  calculated  as  if  the 
interaction  of  the  first  molecule  with  the  others  were  multi¬ 
plied  by  a  coupling  parameter  between  0  and  1.  This  average 
energy  is  in  fact  the  derivative  of  the  cluster  free  energy 
with  respect  to  the  coupling  parameter.  Integration  of  the 
average  energy  with  respect  to  coupling  parameter  then  gives 
the  difference  in  free  energy  of  clusters  of  n  -  1  molecules 
and  clusters  of  n  molecules  (when  an  additional  correction 
is  made  for  the  difference  in  the  size  of  the  containment 
sphere  used  in  the  definition  of  n  molecule  clusters  and  n  -  1 
molecule  clusters). 

The  first  method  depends  for  its  accuracy  on  the  validity 
of  the  assumptions  used  in  the  analytical  calculation  of  the 
free  energy  of  the  low  temperature  cluster.  The  second  and 
third  methods  rely  solely  on  the  computer  simulation  tech¬ 
niques  and  thus  are  clearly  to  be  preferred.  The  cost  and 
difficulty  of  performing  such  simulations  has  decreased 
dramatically  in  the  last  decade,  and  in  the  future  computer 
simulation  methods  will  be  the  most  practical  and  accurate 
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ways  of  evaluating  cluster  free  energies. 

In  this  paper  we  present  a  new  method  of  evaluating 
the  free  energies  and  equilibrium  constants  of  clusters  of 
atoms  and  molecules.  The  method  is  based  on  evaluating  the 
energy  of  a  cluster  as  a  function  of  temperature  using  a 
variation  of  the  molecular  dynamics  computer  simulation  tech¬ 
nique  and  then  integrating  with  regard  to  temperature.  The 
state  of  infinite  temperature  is  used  to  evaluate  the  con¬ 
stant  of  integration  analytically  without  any  approximations. 
Provided  that  the  intermolecular  or  interatomic  potential 
satisfies  certain  conditions,  namely  that  the  potentials  be 

finite  at  all  nonzero  distances  and  diverge  no  more  strongly 
-3 

than  r  at  short  distances,  the  integrand  in  the  tempera¬ 
ture  integration  does  not  diverge  as  the  temperature 
approaches  infinity,  and  thus  the  integration  to  the  infinite 
temperature  state  can  easily  be  performed.  The  method  is 
applied  to  the  calculation  of  the  equilibrium  constants  for 
the  formation  of  clusters  containing  two  to  five  water  mole¬ 
cules  . 

In  Section  II,  we  discuss  the  theory  of  the  method, 
using  Hill’s  physical  cluster  theory  as  a  starting  point. 

In  Section  III  we  discuss  the  theory  of  the  computer  simula¬ 
tion  algorithm  we  use.  Section  IV  discusses  the  intermolecular 
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and  intramolecular  potentials  for  water  that  we  used. 
Section  V  discusses  the  methods  we  used  to  analyze  the 
results  and  verify  their  correctness.  Section  VI  contains 
the  results  for  clusters  of  two  to  five  water  molecules. 
Section  VII  contains  a  brief  discussion  of  the  method. 


II.  Formal  Theory  of  Equilibrium  Constants 


for  Cluster  Formation 


A.  Hill's  general  theory 

Hill's  physical  cluster  theory^  provides  the  statistical 
mechanical  basis  for  calculating  equilibrium  constants  for 
cluster  formation.  His  theory  provides  no  unique  choice  of 
definitions  for  clusters.  Instead  it  provides  a  set  of  easily 
satisfied  conditions  that  the  cluster  definitions  must  satisfy, 
and  for  any  such  set  of  conditions  it  gives  precise  expres¬ 
sions  for  the  cluster  equilibrium  constants.  We  will  apply 
his  general  formulation  to  the  problem  of  water  clusters. 

We  use  classical  mechanics  to  describe  the  nuclear  motions 
in  the  water  molecule.  Let  denote  a  complete  set  of  posi¬ 
tion  and  momentum  variables  for  a  single  water  molecule.  The 
canonical  partition  functions  for  N  water  molecules  in  volume 
V  is 

Q  (P,V)  =  (N'.2Nh^N)  1  J*  dxNexp[ -pH(xN)] .  (2.1) 

v  V 

The  factors  of  two  are  the  symmetry  numbers  for  the  molecules. 
To  apply  Hill's  theory,  one  must  find  a  way  to  decompose  QN 
so  that: 


QN 


(2.2) 
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Here  N  is  a  set  of  nonnegative  integers  N^,  N^,  .  .  . ,  . 

refers  to  the  number  of  clusters  of  i  molecules.  These 
numbers  must  satisfy  the  condition 

ilSNs  -  N-  <2-3> 

The  prime  in  the  sum  in  Eq.  (2.2)  denotes  that  the  sum  is 
over  all  sets  of  nonnegative  values  of  N  that  satisfy  condi¬ 
tion  (2.3).  There  are  many  ways  to  make  such  a  decomposition. 
The  physical  meaning  of  QN  is  that  it  represents  the  canonical 
partition  function  of  a  system  that  contains  N  molecules  that 

exist  as  N.  monomers,  N_  dimers,  ...  N  clusters  of  s  mole- 
1  2.  s 

cules ,  . . . ,  etc . 

The  partition  functions  for  systems  that  contain  only 
one  cluster  play  an  important  role  in  the  theory.  Departing 
from  Hill's  notation,  let  us  define 

q(1)  -(W.((5-v) 

««<2) 

0(n)  -'»oo...i..:te’v>*  "=3  <2-4) 

is  the  canonical  partition  function  for  one  cluster  of 
n  molecules.  According  to  most  reasonable  definitions  of 
clusters,  n  molecules  will  be  in  a  cluster  only  if  they  are 
separated  by  microscopic  distances  that  are  independent  of 
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the  size  of  the  volume  V.  This  leads  to  being  of 

order  V  for  large  V. 

Let  denote  the  average  number  density  of  clusters 
of  n  molecules.  Hill's  theory  provides  a  formula  for  pR 
that  is  an  expansion  in  powers  of  the  density.  To  lowest 
order  in  density,  his  result  is 

*  Kn(P) 

where 

Kn(p)  =  (Q(n)/V)/(Q(1)/V)n.  (2.5) 

The  important  assumptions  made  in  deriving  this  result  are 
that  the  partition  function  can  be  decomposed  according  to 
Eq.  (2.2)  for  all  N  and  that  Q^nVv  approaches  a  limit  as 
V  -  ® . 


B.  One  realization  of  Hill's  theory 

One  way  of  accomplishing  the  decomposition  in  Eq.  (2.2) 
is  to  decompose  phase  space  for  each  value  of  N  into  dis¬ 
joint  regions  such  that  in  each  region  the  states  have  a 
well  defined  set  of  values  of  the  cluster  numbers  N.  In 
particular,  if  we  can  define  the  set  of  functions  C^(x^) 
for  all  N  and  N  that  satisfy  Eq.  (2.3)  such  that 
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(2.6) 


N 

1  if  the  phase  point  x  has  cluster 
numbers  N 

0  otherwise 


and  such  that 

rN  CN(£N)  =  (2.7) 

for  all  N  and  then  Eq.  (2.2)  is  satisfied  if  we  let 

Q  =  (N'^W1  f  dxNC  (xN)exp[-3H(xN)].  (2.8) 

=-  y  £: 

To  define  these  functions  it  is  necessary  to  have  a  pro¬ 
cedure  for  assigning  a  unique  set  of  cluster  numbers  N  to 

N 

every  phase  point  x  for  an  N  molecule  system. 

Let  us  write 


C(1)(xX)  =  1 


C(2)(x2)  =  CQ1(x2) 


(2.9) 


C  ^ (x  )  ~  CQ0  #  # i (x  ) ,  n  a  3 . 


(2.10) 


The  procedure  we  use  for  defining  the  functions  in 
(2.6)  is  the  following.  We  imagine  that  a  necessary  (but  not 
sufficient)  condition  for  n  molecules  to  form  a  cluster  of 
n  molecules  is  that  the  center  of  mass  of  each  of  the  mole¬ 
cules  be  less  than  a  distance  R  from  their  mutual  center  of 


mass.  We  decide  on  the  set  of  values  of  R  for  all  n.  We 

n 

use  the  following  algorithm  for  deciding,  for  a  phase  point 
N 

x  ,  how  many  clusters  of  each  size  are  present. 

1.  Let  n  =  N. 

2.  If  n  =  1,  let  the  molecule  be  regarded  as  a 
cluster  of  1  and  stop. 

3.  If  n  >  1,  see  if  there  are  any  sets  of  n  molecules 

that  are  all  within  a  distance  R  of  their  mutual 

n 

center  of  mass . 

4.  If  there  is  no  such  set,  let  n  =  n-1  and  go  to  2. 

5.  If  there  is  one  such  set,  call  these  molecules  a 
cluster  of  n,  remove  them  from  further  considera¬ 
tion,  let  n  «  n  -  1  and  go  to  2. 

6.  If  there  is  more  than  one  such  set,  find  the  set 
of  n  that  is  most  compact,  in  the  sense  that  the 
largest  distance  of  a  molecule  from  the  mutual 
center  of  mass  is  smaller  for  that  set  than  for 
any  other  set.  Call  this  most  compact  set  a 
cluster  of  n,  remove  them  from  further  considera¬ 
tion,  keep  n  unchanged  and  go  to  3. 

Except  for  a  set  of  phase  points  of  measure  zero  (those  points 
where  there  are  groups  of  the  same  n  that  are  equally  compact), 
this  algorithm  gives  a  precise  way  of  defining  how  many  clusters 
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of  each  type  exist  for  any  phase  point,  and  thus  it  provides 
a  way  of  partitioning  N  molecule  phase  space  in  a  way  needed 
for  the  validity  of  Eqns.  (2.6)  and  (2.7). 

It  follows  from  this  procedure  that 

C^n^(xn)  =1  if  the  center  of  mass  of  each  mole¬ 
cule  is  less  than  R  from  the  center 

n 

of  mass  of  the  set  of  n  molecules 
=  0  otherwise.  (2.11) 

Combining  Eqs.  (2.4),  (2.5),  (2.8),  (2.9),  (2.10),  and 
(2.11),  we  obtain 

Kn(P}  -  In0)/[I10)]n  (2.12) 

where 

I  (P)  =  (n!)  1  J*  dx11  (xn)exp[-pH(xn)] .  (2.13) 

n  y 

C.  Temperature  dependence  of  the  equilibrium  constants 
The  integrals  In(P)  cannot  be  evaluated  directly  by 
computer  simulation,  but  their  logarithmic  derivative  with 
respect  to  p  can  be  evaluated.  From  Eq.  (2.13)  it  follows 
that 

d  In  In/dp  -  -<H(xn)>N 
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where  the  average  on  the  right  side  is  defined  by 


<f(xn)>  =  J  dxnf (xn)exp[-pH(xn) ]C^n^ (xn) 

V 

^  [J  dxnexp[-(3H (xn)]C^n^ (xn)J  1 


It  follows  from  Eq.  (2.12)  that 

d  In  Kn(e)/d£  =  -<H(xn))n  +  n<H(x1)>1.  (2.14) 

The  Hamiltonian  H  is  a  sum  of  kinetic  energy  T  and  potential 
energy  U, 

H(xn)  -  T(xn)  +  U(xn) 


where  U  includes  both  intermolecular  and  intramolecular 
potential  energy.  The  contributions  of  T  to  the  two  terms  on 
the  right  of  (2.14)  cancel  each  other.  If  we  define 

Un(p)  -  <U(xn)>n  (2-15) 

we  have  the  following  result  for  the  temperature  derivative 
of  In  K 

n 


d  In  Kn(p)/d3  -  -Un(3)  +  nUxO). 


Hence, 


In  Kn(p) 


P 

In  K  (0)  -  J  d£J'[U  (£')  -  nU.Oi')] 
n  ©  n 


(2.16) 
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The  first  term  on  the  right  can  be  evaluated  analytically, 
and  the  integrand  on  the  right  can  be  evaluated  by  performing 
molecular  dynamics  simulations.  This  is  the  basic  expression 
we  use  for  the  evaluation  of  the  cluster  formation  equilib¬ 
rium  constants. 

To  evaluate  the  high  temperature  limiting  behavior  of 
In  Kn>  we  note  that  at  high  temperature  the  Boltzmann  factor 
can  be  replaced  by  unity  and  we  have 

Kn(0)  =  InC°)/[I1(°)3n 

■=  (n'.V)"1  J*  dxnC(n)(xn)/[v"1  J*  dx 
V  ~  V 

For  each  molecule  we  convert  the  variables  of  integration  to 
the  momentum  variables,  a  center  of  mass  variable,  and  a 
set  of  internal  coordinates.  The  function  depends  only 

on  the  center  of  mass  variables.  The  momentum  and  internal 
coordinate  integrations  factor  out  of  the  numerator  and 
denominator  and  cancel  each  other,  leaving 

K  (0)  -  (n'.V)”1  J  drnC(n)(rn) 
n  V 

where  the  rn  are  the  center  of  mass  variables.  This  integral 

9 

has  been  evaluated  by  Lee  e_t  al..  The  result  is 
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(2.17) 


Kn(0)  =  a(n)n3/2[A7rR^]n"1/n! 

where  a(n)  is  a  function  of  n  whose  values  range  from  2.437 
to  2.974.  For  n  =  2  to  5,  the  values  are  2.828,  2.436,  2.590, 
and  2.661,  respectively.  (The  values  reported  in  Table  I  of 
Reference  9  are  incorrect.) 

A  potential  complicating  feature  of  the  beta  integra¬ 
tion  in  Eq.  (2.16)  is  that  for  some  inter-molecular  potentials 
the  integrand  diverges  as  3  -  0.  (The  divergence  is  integra- 
ble,  however.)  If  the  integrand  is  very  large,  it  is  diffi¬ 
cult  to  evaluate  accurately  by  computer  simulations,  since  in 
general  the  statistical  noise  in  the  answer  will  also  be  large. 
This  divergence  exists  for  potentials  that  increase  as  r“(3+n) 
as  the  interatomic  or  intermolecular  distance  r  decreases  to 
zero,  for  n  £  0.  For  the  Lennard- Jones  potential,  for  example, 
the . integrand  would  diverge.  The  potentials  we  use  for  water 
diverge  only  as  1/r  for  small  r,  and  the  integrand  approaches 
a  finite  value  as  p  -  0. 
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D.  Choice  of  the  and  the  relationship  between  cluster 

formation  constants  and  virial  coefficients 

According  to  the  definition  of  cluster  that  we  use,  a 

necessary  condition  for  n  molecules  being  a  cluster  of  n 

molecules  is  that  they  all  lie  within  a  sphere  of  radius 

centered  at  their  mutual  center  of  mass.  We  will  refer  to 

this  sphere  as  the  "cluster  containment  sphere."  The  choice 

of  the  containment  sphere  radius  for  each  value  of  n  is 

arbitrary;  for  any  set  of  choices.  Hill's  formal  theory 

provides  expressions  for  the  equilibrium  constants  for  cluster 

formation.  When  the  various  R^  satisfy  certain  conditions, 

it  can  be  shown  that  the  equilibrium  constants  for  dimer  and 

trimer  formation  are  simply  related  to  the  second  and  third 

virial  coefficients.  Here  we  discuss  our  choice  of  the  R 

n 

and  the  relationship  of  equilibrium  constants  and  virial  coef 
ficients . 

For  most  of  our  calculations,  we  chose  R  using  the 

q 

criterion  of  Lee  e_t  al.  i.e.  we  required  that  the  volume  of 
the  sphere  for  clusters  of  n  molecules  be  five  times  the 
volume  occupied  by  n  molecules  in  the  liquid  at  standard 
temperature  and  pressure.  Thus 

R  -  3.292  n1/3  k. 
n 
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The  factor  of  five  is  admittedly  arbitrary,  but  Lee  e£  al . 
demonstrated  that  the  free  energy  of  clusters  of  atoms  is 
independent  of  the  volume  of  the  containment  sphere  for  low 
enough  temperatures.  For  clusters  of  three  water  molecules, 
we  have  also  performed  calculations  with  a  different  choice  of 
and  have  verified  that  the  equilibrium  constants  are 
insensitive  to  the  volume  of  the  containment  sphere  at  low 
temperatures . 

The  equilibrium  constant  for  dimer  formation  can  be 
related  to  the  second  virial  coefficient  for  the  gas  under 
certain  conditions.  When  two  molecules  are  in  a  cluster, 
they  are  less  than  a  distance  2R£  apart.  Thus  the  dimer  con¬ 
stant  is  not  affected  by  the  form  of  the  intermolecular  poten¬ 
tial  for  separations  larger  than  this.  If  the  potential  is 
zero  at  such  large  distances,  it  is  straightforward  to  show, 
using  Eqns.  (2.1)  and  (2.8)  that 

Q2  -  Q(2)  +  [V- y(2R2)3]Q2/2V. 

Using  the  standard  relationship  between  the  second  virial 

coefficient  and  the  partition  functions  for  one  and  two 
11 


molecules 


we  find 


B2  =  16ttR^/3  -  Kr  (2.18) 

Thus,  if  K2  is  calculated  for  a  particular  choice  of  R2  and 
a  particular  intermolecular  potential,  this  equation  shows 
that  the  result  can  be  used  to  obtain  the  second  virial  coef¬ 
ficient  for  a  gas  whose  intermolecular  potential  is  zero  for 
separations  greater  than  2R^  and  equal  to  the  potential  used 
in  the  calculation  for  separations  less  than  2R2- 

A  similarly  simple  relationship  for  the  third  virial 
coefficient  cannot  be  obtained  for  an  arbitrary  choice  of  R^. 
However,  if 

R3  -  2R2  (2.19) 

and  if  the  range  of  the  potential  is  equal  to  or  less  than 
2R2,  we  can  obtain  a  simple  result.  (A  simple  result  can 
also  be  obtained  if  R^  >  2R2>  but  we  shall  not  discuss  this 
case.)  When  these  restrictions  are  satisfied,  if  three 
molecules  are  located  so  that  there  are  two  (or  three)  inter¬ 
molecular  distances  less  than  2R2>  the  three  molecules  must 
form  a  trimer.  The  points  in  three  molecule  configuration 
space  can  then  be  divided  into  the  following  sets:  A.  points 
in  which  the  three  molecules  form  a  trimer;  B.  points  in 


which  the  three  molecules  do  not  form  a  trimer  but  two 


molecules  form  a  dimer;  and  C.  points  in  which  the  three 
molecules  do  not  form  any  dimer  and  no  intermolecular  dis¬ 
tance  is  less  than  2R2>  In  set  C,  all  the  intermolecular 
interactions  are  zero.  In  set  B,  the  molecule  that  is  not 
in  the  dimer  must  be  separated  from  the  center  of  mass  of  the 
other  two  by  at  least  (3/2)R^>  which  is  equal  to  3R2 >  to  avoid 
formation  of  a  trimer.  Since  the  other  two  are  separated  by 
a  distance  of  less  than  2R2,  this  latter  condition  also 
guarantees  that  the  third  molecule  is  further  than  R2  from 
each  of  the  others  and  hence  does  not  interact  with  the  other 
two.  It  is  then  straightforward  to  show  that 

Q3  -  Q(3)  +  [V-  (47r/3)(3R2)3]Q1Q(2)/V  +  aQ^ 

where  the  three  terms  on  the  right  correspond  to  the  regions 
A,  B,  and  C,  respectively.  The  quantity  a  is  an  integral 
that  depends  on  V  and  on  R2  but  is  independent  of  p.  Using 
the  standard  relationship  between  the  third  virial  coefficient 
and  the  two  and  three  molecule  partition  functions,^ 

B3  -  -V2[2(Q3/Q^) -4(q2/q^)  +  2(Q2/Q^) -y], 
and  Eq.  (2.18),  we  find 

B3  -  -2K3'  +  4B2  -  72ttR2B2  + 


{(j  -  2a)v2  -  32*^73  +  384ir2R2} 


where  the  prime  in  the  trimer  equilibrium  constant  denotes 
that  it  is  defined  using  the  nonstandard  trimer  containment 
sphere  radius  given  in  Eq.  (2.19).  The  integral  a  can  be 
evaluated  tediously.  A  simpler  procedure  is  to  use  the  fact 
that  for  p  =  0  the  virial  coefficients  are  zero  and  the  equi¬ 
librium  constants  are  given  by  Eq.  (2.17).  Then  we  obtain 

B3  *=  -2K.J  +  4B2  -  72ttR^B2  +  480tt2R*.  (2.20) 

This  relationship  between  the  second  and  third  virial  coef¬ 
ficient  and  the  dimer  and  trimer  equilibrium  constants  is 
valid  when  Eq.  (2.19)  holds  and  when  the  range  of  the  inter- 
molecular  potential  is  equal  to  or  less  than  2R2» 


III.  Theory  of  the  Simulation  Algorithm 


A.  Description  and  justification  of  the  algorithm 

To  evaluate  Un((3),  defined  in  Eq.  (2.15),  we  must 
evaluate  the  average  of  the  n  molecule  potential  energy  over 
a  distribution  of  n  molecule  phase  points,  where  the  distri¬ 
bution  function  in  n  molecule  phase  space  is  proportional  to 

exp[-pH(xn)]C^(xn) .  (3.1) 

The  function  can  be  regarded  as  the  Boltzmann  factor 

for  a  hard  potential  that  keeps  the  n  molecules  within  the 
region  of  space  in  which  they  form  an  n  molecule  cluster. 

If  we  define 

W(n)(xn)  -  0  ifC(n)(xn)-l 

-  •  if  C(n)(xn)  -  0, 

then  the  distribution  function  in  (3.1)  can  be  written  as 

exp[-p{H(xn)  +  (x11)}]  ,  (3.2) 

which  is  a  Boltzmann  factor  for  the  combination  of  the 
molecular  Hamiltonian  and  the  wall  potential.  From  (2.11) 
it  can  be  seen  that  the  wall  potential  is  a  function  of  the 
center  of  masses  of  the  n  molecules  and  is  independent  of 


r 


their  internal  coordinates.  Since  the  wall  potential  depends 

only  on  the  relative  positions  of  the  centers  of  mass,  it  is 

a  momentum  conserving  interaction.  When  the  molecules  are 

within  the  cluster  containment  sphere  of  radius  R  centered 

n 

at  their  mutual  center  of  mass,  the  potential  is  zero,  but  it 
becomes  nonzero  and  infinite  when  the  center  of  mass  of  any 
of  the  molecules  touches  the  surface  of  the  cluster  contain¬ 
ment  sphere.  When  this  happens,  the  relative  velocity  of  the 
center  of  mass  of  that  molecule  with  respect  to  the  mutual 
center  of  mass  of  the  remaining  n  -  1  molecules  is  reversed  in 
direction  and  unchanged  in  magnitude.  In  effect,  the  molecule 
that  hits  the  surface  of  the  containment  sphere  is  specularly 
reflected  from  the  surface  and  the  velocity  of  each  of  the 
other  molecules  is  changed  by  the  amount  required  to  conserve 
total  momentum.  Each  of  the  other  molecules  suffers  the  same 
change  in  center  of  mass  velocity,  and  the  positions,  internal 
coordinates,  and  internal  velocities  of  all  the  molecules  are 
unaffected.  The  collisions  with  the  wall  conserve  kinetic 
energy. 

To  generate  a  set  of  phase  points  for  an  n  molecule 
cluster  that  are  distributed  according  to  the  Boltzmann  factor 
in  (3.2),  we  use  the  following  procedure,  which  is  a  combina¬ 
tion  of  the  molecular  dynamics  and  Monte  Carlo  simulation 
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methods.  We  calculate  a  molecular  dynamics  trajectory  for 
a  system  of  n  molecules  subject  to  the  Hamiltonian  H  and  to 
the  wall  potential  W.  At  regular  intervals  of  time,  we  replace 
the  momenta  of  each  of  the  atoms  by  a  momentum  chosen  at  random 
from  the  Boltzmann  distribution  appropriate  for  the  mass  of 
the  atom  and  for  the  temperature  of  interest.  This  latter  pro¬ 
cedure  can  be  regarded  as  subjecting  the  atoms  to  stochastic 
collisions  with  a  heat  bath  that  has  the  temperature  of  inter¬ 
est.  Using  the  theory  of  Markov  processes,  it  is  possible  to 
show  that  this  combination  of  Hamiltonian  dynamics  and  sto¬ 
chastic  collisions  is  a  Markov  process  and  that  a  trajectory 
for  the  process  has  the  property  that  the  time  average  of  a 
mechanical  property  over  a  trajectory  is  equal  to  the  statis¬ 
tical  average  of  the  same  property  over  a  distribution  function 
proportional  to  (3.2),  provided  the  trajectory  is  long  enough. 
(The  details  of  the  proof  are  similar  to  those  of  a  proof  in 
Reference  12  and  they  will  be  omitted  here.) 

Three  important  assumptions  must  be  made  to  apply  the 
theorem.  First  it  must  be  assumed  that  under  the  motion 
generated  by  the  Markov  process  a  system  can  in  principle  get 
from  any  phase  point  to  any  other  phase  point  in  a  finite 
amount  of  time.  (In  other  words,  the  motion  must  be  ergodic.) 
The  second  assumption  is  that  the  trajectory  used  in  the 
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calculation  is  long  enough  so  that  the  trajectory  actually 
samples  all  the  important  parts  of  phase  space.  The  third 
assumption  is  that  the  trajectory  is  long  enough  that  the 
statistical  noise  in  the  calculation  is  averaged  away. 

For  dense  liquids  it  is  commonly  believed  that  the  first 
assumption  is  correct,  and  it  is  very  reasonable  to  expect 
that  it  is  correct  for  this  type  of  simulation  of  a  cluster. 
The  second  assumption  can  be  tested  by  starting  two  trajec¬ 
tories  in  very  different  parts  of  phase  space  and  seeing  if 
they  both  give  the  same  time  average.  The  statistical  error 
in  a  calculation  can  be  estimated  on  the  basis  of  physical 
reasoning  or  from  the  results  of  the  simulation.  In  the  next 
part  of  this  section  we  discuss  how  the  statistical  error  can 
be  estimated. 

B.  Estimation  of  statistical  error 

Let  A  represent  a  dynamical  variable.  Ensemble  averages 
of  any  quantity  over  the  correct  distribution  will  be  denoted 
by  angular  brackets.  Time  averages  over  a  particular  trajec¬ 
tory  will  be  denoted  by  overbars.  The  time  average  of  A  is 

1  T 

A  =  T  1  J  dt  A(t) , 

0 

where  A(t)  is  the  value  of  A  at  time  t  on  the  particular  tra¬ 
jectory.  The  theorem  mentioned  above  states  that  under  the 
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appropriate  conditions 


A  -  <  A)  as  t  ~  ®  . 

Let 

A (t)  =  A(t)  -  <  A>  . 

Then  the  error  made  in  assuming  A  and  (A)  are  equal  is 

-1  T 

A  -  <A>  =  T  J  dt  Lit). 

0 

The  square  of  the  error  is 

T  T 

(A  -  <  A)  ) 2  -  T-2  J  dt.  |  dt.  A(t.)A(t  ). 

o  x  o  c  * 

We  can  estimate  the  square  of  the  error  by  calculating  the 
ensemble  average  of  both  sides  of  this  equation;  i.e.  we 
calculate  the  average  over  an  ensemble  of  trajectories  start¬ 
ing  at  all  points  in  phase  space  and  appropriately  weighted 
by  the  probability  distribution  of  the  initial  state.  We 
obtain 

T  T 

<(A-<A>)2)  «  T"2  J  dt,  J  dt  <A(t,)A(t2». 

0  0  z 

The  average  in  the  integrand  is  the  ensemble  averaged  two 

time  correlation  function  of  the  fluctuation  of  A  from  its 
ensemble  average.  For  the  ensemble,  such  a  correlation 
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I 


function  depends  only  on  the  time  interval  l^-tgl,  and  we 
expect  the  correlation  function  to  vanish  as  1 1.  -  t  I  -  ® 
Thus  we  can  write 

<A(t1>A(t2)>  -  <A2>g(t1~t2) 

where 

g(0)  =  1 

g(t)  -  0  as  t  -  ±  ®. 


We  obtain 

<(A-<A>)2)  =  2<A2>tc/T 


(3.3) 


where 


T 

C 


I  dt  g(t) 
0 


In  obtaining  Eq.  (3.3)  we  have  assumed  that  T  is  large  com¬ 
pared  with  the  times  for  which  g(t)  is  nonzero.  The  quantity 
tc  can  be  regarded  as  the  correlation  time  for  fluctuation 
of  A  about  its  ensemble  average.  (If  g(t)  were  an  exponential 
function  of  t,  tc  would  be  the  time  constant  for  the  exponen¬ 
tial.)  Let  us  define 


N 


c 


N  is  the  number  of  correlation  times  contained  within  the 
c 

duration,  T,  of  the  trajectory.  Hence 
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(3.4) 


<  (A- (A>)2)1/2  =  <aV/2(2/Nc)1/2. 

This  is  an  estimate  of  the  root  mean  square  error  made  by 
assuming  that  A  obtained  from  one  trajectory  is  equal  to  (A). 

To  estimate  the  right  side  of  this  equation  it  is  reason 
able  to  use  the  results  of  the  one  trajectory  and  assume 

<A2)  =  <  (A  -  ( A)  )2>  «  (A-S)2  (3.5) 


and 

<A2>g(t)  -  < A (0)a (t)> 

_1  T-t 

*  (T-t)  J  dt  [A(t  )-A][A(t.+t)-A]  (3.6) 

0 

i.e.  to  assume  that  the  trajectory  averaged  correlation  func¬ 
tion  of  fluctuations  of  A  from  its  trajectory  average  is 

approximately  equal  to  the  ensemble  average. 

2 

The  time  average  (A  -  5)  is  easy  to  evaluate  by  perform¬ 
ing  the  indicated  integration  of  fluctuations  of  A.  (Alterna 
tively,  when  A  is  a  potential  energy,  as  it  is  in  our  case, 
the  ensemble  average  squared  fluctuation  in  the  potential 
energy  is  simply  related  to  the  temperature  derivative  of 
the  average  potential  energy,  which  can  be  estimated  from 
numerical  differentiation  of  the  averages  with  regard  to 
temperature.)  The  correlation  function  on  the  right  side  of 


30 


the  last  part  of  (3.6)  is  more  difficult  to  evaluate,  but  it 
can  be  obtained  from  the  results  of  the  simulation.  Alterna¬ 
tively  one  could  use  physical  reasoning  to  estimate  t  . 

The  quantity  can  be  interpreted  as  the  time  over  which 
fluctuations  in  A  are  correlated  rather  than  uncorrelated.  If 
A  represents  intermolecular  potential  energy,  as  in  the  problem 
at  hand,  then  there  are  some  parts  of  phase  space  where  A  is 
large  and  negative  (namely  in  hydrogen  bonding  configurations), 
some  parts  where  A  is  small  (namely  when  the  molecules  are  far 
apart),  and  other  parts  where  A  is  large  and  positive.  At  any 
temperature,  a  certain  range  of  energies  is  likely  to  be  import¬ 
ant  in  the  equilibrium  distribution.  Then  at  any  temperature, 
tc  should  be  approximately  the  time  scale  for  moving  from  the 
low  energy  to  the  high  energy  parts  of  phase  space  that  are 
important  at  that  temperature . 

C.  Choice  of  interval  between  stochastic  collisions 

An  important  parameter  in  the  calculation  is  the  time 
interval  between  the  stochastic  collisions  suffered  by  the 
molecules.  For  infinitely  long  trajectories,  the  value  of 
the  time  interval  is  irrelevant:  correct  average  potential 
energies  will  be  obtained  for  any  choice  of  interval.  For 
trajectories  of  finite  length,  the  value  of  the  interval  has 
an  important  effect  on  the  statistical  error  of  the  calculation. 
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and  thus  it  must  be  chosen  carefully. 

The  stochastic  collisions  equilibrate  the  kinetic 

energy  and  the  total  energy  of  the  water  molecules  to  the 

values  typical  of  the  temperature  of  interest.  Taking  only 

this  into  account,  one  might  be  tempted  to  make  the  collision 

frequency  very  large  to  insure  rapid  equilibration.  However, 

this  would  impede  the  motion  of  the  molecules  in  configuration 

space.  In  effect,  all  motion  would  be  diffusion  limited.  As 

discussed  above,  it  is  important  that  the  molecules  be  able 

to  visit  all  the  appropriate  high  energy  and  low  energy  parts 

of  configuration  space  in  order  to  obtain  good  statistical 

averages  from  the  trajectory.  If  the  stochastic  collisions 

are  too  frequent,  Tc  will  be  too  large,  thereby  making  the 

statistical  error  large.  If  we  can  estimate  tc  in  the  absence 

of  stochastic  collisions,  it  is  worthwhile  to  choose  the  time 

interval  between  stochastic  collisions  to  be  no  shorter  than 

this  estimate  of  t  .  Stochastic  collisions  that  are  this 

c 

infrequent  will  not  inhibit  the  motion  of  the  system  from  high 
to  low  energy  parts  of  configuration  space.  Thus,  in  our 
calculations  on  a  particular  system  at  a  particular  tempera¬ 
ture,  we  make  an  estimate  (on  the  basis  of  physical  reasoning 
or  previous  calculations)  of  what  t will  be  for  that  system 
and  we  choose  the  interval  between  stochastic  collisions  to 


be  that  estimate. 


IV.  The  Molecular  Interactions 

For  this  simulation  of  water  clusters  the  intramolecular 
potential  was  assumed  to  be  a  sum  of  atom-atom  potentials. 

For  the  two  intramolecular  atom-atom  potentials,  we  used 
harmonic  potentials  with  equilibrium  bond  lengths  and  spring 
constants  which  correspond  to  the  minima  and  curvatures  of 
the  improved  central  force  potentials  of  Stillinger  and 
Rahman^: 

V0H(r)  "  f  0-147. 6)  (r  -  0.9584  Jt)2  kcal/mole 

VHH(r)  =  f<257.3)(r  -  1.5151  A)2  kcal/mole. 

The  minima  of  these  two  potentials  are  those  which  reproduce 
the  correct  geometry  of  an  isolated  water  molecule.  The 
curvatures  were  selected  to  reproduce  the  asymmetric  stretch 
frequency  in  D2O  and  to  equalize  the  fractional  errors  in 
the  frequencies  of  the  other  two  vibrational  modes  of  the 
same  molecule . 

To  define  the  intermolecular  potential,  let  r^  repre¬ 
sent  the  position  of  the  ith  atom  on  the  jth  molecule.  The 
first  atom  on  each  molecule  is  the  oxygen  atom  and  the  others 
are  hydrogen  atoms.  We  define  a  set  of  intermolecular  atom- 
atom  interactions  V^Cr),  V^Cr),  and  V^Cr),  and  a  switching 
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function  S(r).  Then  the  intermolecular  interactions  are  of 
the  form: 


U  ^11  >~21  *~31  *~12  *~22  ’~32  * 

=  S(ten  “£i2I>.  ^  Vij(teii 

i,j=l  J  _  J 


where  the  subscripts  on  the  V„  functions  should  be  00,  OH, 

or  HH,  as  is  appropriate  to  the  nature  of  the  atoms  i  and  j . 

For  these  atom-atom  intermolecular  interactions,  we  chose  to 

14 

use  those  of  Watts  which  were  optimized  to  reproduce  the 
temperature  dependence  of  the  second  virial  coefficient  of 
water  in  the  temperature  range  for  which  experimental  data 
was  available.  The  switch  function,  S(r),  is  designed  to 
switch  the  Watts  potential  off  smoothly  as  a  pair  of  molecules 
separates.  The  switch  function  has  been  used  before^  with 
excellent  results.  The  form  of  the  S(r)  is: 


S(r) 


for  0  <  r  <  rT 


s(r  )  for  rL  <  r  <  ry 
0  for  r^  <  r. 


2  2 
where  s(r  )  is  a  fifth-order  polynomial  function  in  r  which 

2 

is  designed  to  have  zero  first  and  second  derivatives  at  r^ 

0% 

and  r|j,  and  which  lets  S(r)  be  continuous.  (We  chose  6.0  k 
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and  5.5  A  for  and  r^»  respectively.)  The  use  of  switch 
functions  is  preferable  to  a  simple  truncation  (or  trunca¬ 
tion  and  shift)  of  the  atom-atom  potentials  especially  when 
the  values  of  these  potentials  or  the  forces  are  large 
where  the  truncation  is  to  occur.  This  is  the  case  for  the 
water-water  potentials  used  here  because  the  three  atom- 
atom  components  are  dominated  by  a  long-range  Coulombic 
term  at  large  distances  even  though  their  sum,  the  full  Watts 
potential,  is  dipolar-dipolar  at  these  distances.  The  use 
of  a  continuous  switch  function  in  the  potential  is  also 
preferable  to  truncation  of  the  force  for  intermolecular  dis¬ 
tances  beyond  a  certain  cutoff  distance.  Truncation  of  the 
force  leads  to  a  force  that  is  not  the  derivative  of  a  poten¬ 
tial  and  hence  leads  to  lack  of  conservation  of  energy  and  a 
secular  heating  of  the  sample . 

The  choice  of  a  purely  harmonic  intramolecular  potential 
was  motivated  by  a  desire  for  simplicity  and  by  the  expecta¬ 
tion  that  cluster  equilibrium  constants  are  more  sensitive  to 
intermolecular  interactions  than  to  intramolecular  inter¬ 
actions.  The  Watts  intermolecular  potential  was  chosen  for 
these  cluster  studies  because  it  has  the  correct  second  virial 
coefficient.  The  switching  off  of  the  Watts  potential  at 
long  distances  is  expected  to  have  a  small  effect  on  the 
properties  of  the  small  clusters  we  are  concerned  with  here. 
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V.  Computational  Procedure,  Error  Analysis 
and  Consistency  Checks 

Clusters  of  one  to  five  molecules  were  simulated  by 

molecular  dynamics  using  the  potentials  described  in  the 

previous  section.  For  each  cluster  size,  the  system  was 

equilibrated  to  several  different  temperatures  through  the 

administration  of  stochastic  collisions  on  every  atom  at 

regular  intervals.  The  potential  energy,  U^,  of  each  n 

molecule  cluster  was  calculated  at  each  molecular  dynamics 

time  step  and  averaged  over  the  length  of  a  run.  Finally, 

Un  -  nU^  was  plotted  and  integrated  as  a  function  of  p 

(=  l/k_T)  to  yield  In  K  (p)  .  This  section  describes  some 
is  n 

aspects  of  the  computational  procedure,  states  the  parameters 
of  the  calculation,  and  discusses  several  consistency  checks 
used  to  verify  the  accuracy  of  the  calculation. 

To  choose  the  time  intervals  between  stochastic  colli¬ 
sions,  we  estimated  the  correlation  time  for  potential  energy 
fluctuations  using  physical  reasoning  and  choose  the  interval 
to  be  that  estimate.  For  monomers,  the  correlation  time  for 
potential  energy  fluctuations  was  estimated  to  be  the  period 

of  the  slowest  vibrational  mode,  the  bending  vibration, 

-14 

2  x  10  sec.  For  clusters,  we  considered  two  possibilities 
for  the  time  scale  for  energy  fluctuations.  The  first  is 
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the  time  for  a  molecule  to  cross  the  containment  sphere.  At 
high  temperatures,  the  molecules  are  not  bound  to  each  other 
and  the  time  scale  for  potential  energy  fluctuations  is 
approximately  the  time  between  collisions  of  a  molecule,  which 
is  approximately  the  sphere  traversal  time.  The  second  is  the 
intermolecular  vibrational  period  for  a  cluster,  which  we 
estimated  as  0.09  psec.  This  is  an  appropriate  estimate  of 
the  correlation  time  at  low  temperatures.  We  actually  used 
the  geometric  mean  of  these  two  estimates  as  the  interval 
between  stochastic  collisions  for  our  simulations  of  clusters 
at  all  temperatures.  Typical  values  of  the  interval  were  0.3 
to  0.6  psec. 

The  time  duration  of  each  simulation  was  chosen  on  the 
basis  of  the  amount  of  statistical  error  that  would  be  tole¬ 
rated.  To  estimate  the  statistical  error  as  a  function  of 
duration,  Eq.  (3.4)  was  used.  The  correlation  time  needed  for 
the  right  side  of  Eq.  (3.4)  was  estimated  to  be  equal  to  the 
geometric  mean  of  the^sphere  traversal  time  and  the  cluster 
breathing  period.  The  mean  square  fluctuation  of  the  poten¬ 
tial  energy  was  estimated  using  Eq.  (5.1)  below  and  an  esti¬ 
mate  of  (27n  -  6)kg/2  as  the  heat  capacity  of  an  n-mer.  (Note 
that  the  statistical  error  estimated  in  this  way  was  done  in 
advance  of  the  simulation  and  the  estimate  was  used  only  for 
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deciding  on  the  length  of  the  simulation.  The  statistical 
errors  reported  in  the  results  section  below  were  obtained 
from  correlation  times  and  mean  square  energy  fluctuations 
obtained  from  the  simulation  data.)  The  durations  varied 
from  five  nanoseconds  (for  the  dimer  at  2000?  K)  to  400  pico¬ 
seconds  (for  the  monomer  at  50°  K) . 

The  statistical  errors  in  the  resulting  potential  energy 
averages  were  calculated  by  computing  the  mean  square  fluc¬ 
tuation  in  the  potential  energy  for  each  simulation  and  the 
correlation  function  of  the  fluctuations  for  some  of  the  simu¬ 
lations  and  using  Eqs .  (3.4)- (3.6).  Typical  values  of  the 
error  were  0.02,  0.2,  0.2,  0.3,  and  0.4  kcal/mole  for 
Uj  (p)  . .  .Uj.  (p) ,  respectively. 

The  choice  of  time  step  for  the  numerical  integration  of 
the  equations  of  motion  is  crucial  for  obtaining  accurate 
results.  The  step  must  be  small  enough  that  the  integration 
algorithm  gives  an  accurate  description  of  the  trajectory,  at 
least  for  the  purpose  of  calculating  potential  energy  averages. 
Too  small  a  value  is  undesirable,  because  for  a  given  amount 
of  computer  time  the  number  of  correlation  times  in  the  tra¬ 
jectory  is  inversely  proportional  to  the  time  step.  We  chose 
the  time  step  in  the  following  way.  We  chose  the  durations 
of  the  runs  so  that  the  statistical  error  in  the  potential 
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energy  averages  would  be  about  0.1  kcal/mole.  We  wanted 
the  time  step  to  be  small  enough  that  the  trajectory  calcu¬ 
lated  for  times  of  the  order  of  the  correlation  time  would 
give  average  potential  energies  that  are  in  error  by  no  more 
than  0.05  kcal/mole.  We  performed  various  test  calculations 
starting  at  the  same  initial  mechanical  state  and  integrat¬ 
ing  the  equations  of  motion  for  the  same  length  of  time 
using  different  time  intervals,  and  we  determined  what  the 
time  step  should  be  to  give  the  desired  accuracy.  As  the 
result  of  these  tests,  we  used  a  timestep  of  0.2  x  10~^  sec 
in  all  calculations  except  those  at  400(f  K,  which  used  0.1 
x  10  sec. 

The  equations  of  motion  were  integrated  using  a  version 
of  the  Verlet  algorithm  discussed  in  the  appendix. 

U^,  the  monomer  energy,  as  a  function  of  T,  was  fit 

3 

very  well  by  a  straight  line  of  slope  —  k  .  This  is  not  sur- 
prising  since  the  atoms  within  a  molecule  interact  by  purely 
harmonic  potentials.  There  was  substantial  deviation  from 
the  straight  line  behavior  only  at  the  one  high  temperature 
of  400CP  K.  (This  may  be  due  to  rotational-vibrational  coup¬ 
ling.)  For  all  lower  temperatures,  the  linear  fit  was  used 
to  obtain  nU^,. 

For  each  cluster  size  the  data  consisted  of  U  -  nU.  at 

n  i- 

12  to  15  values  of  the  temperature  (see  Figure  1).  Since 
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it  is  the  area  under  this  curve  as  a  function  of  £  =  1/k  T 

D 

that  gives  d  In  K^/dp,  the  temperature  values  selected  were 
more  or  less  evenly  spaced  in  £  and  ranged  from  £  *  0.25 
(kcal/mole)-1  (~  2000p  K)  to  £  =  3.0  (kcal/mole)-1  (~  17CfK). 

For  each  cluster  size,  the  function  U  -  nU,  was  fitted 

n  1 

to  an  analytic  form  to  facilitate  the  integration  of  In  Kn(P)* 
The  fitting  procedure  called  for  a  simple  analytic  form,  so 
that  the  integration  would  be  easy  to  perform,  as  well  as 
some  sort  of  smoothing,  to  decrease  the  effect  of  statistical 
noise . 

The  fitting  procedure  chosen  used  a  cubic  spline.  For 
the  m  values  of  and  y^,  the  fit  is  determined  by  the  m-1 
sets  of  the  four  coefficients  that  define  the  cubic  polynomial 
in  each  interval.  The  cubic  spine  fit  is  further  required  to 
be  continuous  and  have  continuous  first  and  second  derivatives 
at  each  x^.  Most  cubic  spline  fits  further  require  the  fit 
to  pass  through  all  the  data  points.  To  allow  for  smoothing, 
this  last  condition  was  relaxed.  Given  a  set  of  error  esti¬ 
mates  6y^  for  each  data  point,  the  procedure  that  was  used^ 
produced  the  unique  cubic  spline,  S(x),  such  that  the  quantity 

Jmdx[S'(x)]2 
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was  minimized  subject  to  the  constraint  that 


m 

I 

i=l 


S(xi) 


2 


<:  P. 


If  P  =  0,  the  cubic  spline  goes  through  the  data  points. 

For  very  large  P,  the  spline  approaches  a  linear  least  squares 
fit. 

The  fitted  function  S(x)  is  unique  for  each  choice  of 
P.  Following  earlier  work,^  P  is  determined  to  be  the  value 
which  allows  the  fit  to  go  "between"  the  data  points  rather 
than  consistently  above  or  below  them.  This  idea  is  made 
more  precise  by  defining  the  residual 

R.  =  y .  -  S(x. ) 

1  l  l 

and  a  correlation 


m-1 

E  R .  R .  . 
i-i  1  1+1 


Q  is  a  measure  of  the  amount  of  smoothing.  For  large  values 
of  P,  adjacent  data  points  tend  to  fall  on  the  same  side  of 
the  fit  and  Q  is  positive.  For  small  values  of  P,  Q  is 
invariable  negative.  Furthermore,  Q  is  usually  a  well-behaved 
function  of  P.  The  value  of  P  that  yields  Q  ■  0  is  taken  as 
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the  value  that  yields  a  fit  that  is  smoothed  enough  to  remove 
statistical  noise  without  removing  the  features  of  interest. 

The  values  used  for  were  the  statistical  error  in 
the  U  .  Note  that  an  acceptable  fit  should  also  satisfy 
<  6y^  for  all  i.  This  condition  was  usually  satisfied. 

Once  the  data  was  plotted  and  the  fit  was  determined, 
it  was  quite  easy  to  integrate  the  d  In  K^/dp  curve  from 
P  =  0  to  any  other  value.  As  explained  in  Section  II,  there 
is  no  divergence  in  -  nU^  near  p  =  0. 

The  remainder  of  this  section  will  discuss  several  con¬ 
sistency  checks  that  were  made  by  comparing  different  cal¬ 
culations  in  order  to  detect  various  types  of  systematic  and 
random  errors . 

One  kind  of  error  possible  in  this  procedure  arises  from 
there  being  two  or  more  low-energy  regions  of  phase  space 
with  different  values  of  the  potential  energy  and  separated  by 
a  large  potential  barrier.  At  low  enough  temperatures,  the 
rate  of  barrier  crossing  will  vanish  and  average  values  of  Ur 
will  be  different  depending  upon  the  side  of  the  barrier  on 
which  the  system  remains  trapped.  (One  possible  example 
might  be  the  existence  of  two  stable  n-molecule  clusters, 
one  in  an  n-membered  ring  and  one  in  an  n-1  membered  ring 
with  a  branch.)  One  way  to  test  for  the  existence  of  such  a 


problem  is  to  do  two  or  more  simulations  from  different 
starting  configurations  and  notice  if  the  results  obtained 
are  consistent  with  each  other--i.e.,  if  each  lies  within 
the  error  ranges  of  the  others.  We  performed  such  tests  at 
low,  medium  and  high  temperatures  for  clusters  of  all  sizes 
considered  to  check  that  the  computed  was  the  same  at  each 
temperature  for  all  starting  configurations.  If  there  were 
low  energy  but  mutually  inaccessible  regions  of  phase  space 
not  sampled  during  our  simulations  they  either  did  not  show 
up  in  our  tests  or  the  values  of  the  energy  in  each  region 
were  very  close.  To  obtain  the  plots  of  Un  -  nU^  which  were 
fitted  and  integrated,  the  values  of  Un  from  different  initial 
conditions  were  averaged  together. 

A  second  consistency  check  is  to  see  if  the  fluctuations 
in  potential  energy  at  a  given  temperature  are  correctly 
related  to  the  heat  capacity  at  that  temperature.  For  a 
canonical  ensemble , 

U2  -  U2  -  -du"- /dp.  (5.1) 

n  n  n 

2 

To  make  this  check  U  was  calculated  and  averaged  at  each 

n 

time  step  so  that  the  left  side  of  Eq.  (5.1)  could  be  deter¬ 
mined.  dU  /dS  was  determined  from  the  cubic  spline  fit  of 
n 

U  -nU,  as  a  function  of  B  and  the  linear  fit  of  U.  as  a  func- 
n  1  1 

tion  of  T.  The  results  of  this  test  for  dimers  and  pentamers 


43 


I 


are  displayed  in  Figure  2,  and  those  for  trimers  and  tetramers 
are  similar.  In  this  figure  the  curves  are  more  accurate  than 
the  dots,  since  each  curve  is  the  smoothed  result  of  calcula¬ 
tions  of  the  average  energy  for  all  simulations  of  an  n  molecule 
cluster,  whereas  each  dot  is  the  result  of  a  calculation  of 
the  energy  fluctuations  for  just  one  simulation.  If  the  dots 
fell  exactly  on  the  curve,  this  would  mean  that  each  simula¬ 
tion  was  long  enough  for  the  mean  squared  fluctuations  in  the 
potential  energy  to  be  calculated  accurately.  The  statistical 
errors  in  average  values  is  usually  smaller  than  the  errors 
in  mean  squared  fluctuations,  and  thus  agreement  of  the  points 
and  the  curve  would  imply  that  the  average  potential  energies 
are  subject  to  very  little  statistical  error.  If  the  points 
fall  near  but  not  exactly  on  the  line,  it  is  still  likely 
that  the  trajectories  are  long  enough  to  give  accurate  averages 
even  though  their  fluctuations  are  in  error. 

At  high  or  low  temperature,  where  the  system  is  behaving 
mostly  as  n  noninteracting  molecules  or  as  an  n  molecule 
cluster,  respectively,  the  agreement  of  the  points  and  the 
curve  is  excellent.  The  agreement  is  much  worse  at  inter¬ 
mediate  temperatures  where  there  is  a  local  maximum  in  the 
cluster  specific  heat.  In  fact,  the  correlation  times  evalu¬ 
ated  from  the  simulations  in  the  intermediate  region  are 
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generally  longer  than  the  prior  estimates  based  on  the 
sphere  traversal  time  and  the  cluster  intermolecular  vibra¬ 
tional  frequency.  This  probably  means  that  the  time  for 
formation  and  breakup  of  clusters  determines  the  correlation 
time.  For  these  temperatures,  the  runs  were  not  long  enough 
to  achieve  the  same  statistical  accuracy  achieved  at  lower 
and  higher  temperatures. 

The  volume  of  the  containment  sphere  should  have  no 
effect  on  either  -  nU^  or  at  low  temperatures,  provided 
the  volume  is  large  enough  not  to  interfere  with  the  motion 
of  the  n  molecule  cluster.  To  test  that  this  was  so,  we  per¬ 
formed  a  new  series  of  simulations  for  the  trimer  using  a 
sphere  that  is  5.33  times  larger  in  volume  than  the  standard 
trimer  containment  sphere.  Although  there  should  be  dif¬ 
ferences  in  and  Kn  at  high  temperature,  these  should  dis¬ 
appear  at  temperatures  low  enough  that  the  cluster  is  compact 
and  the  molecules  stay  away  from  the  wall  of  the  containment 
sphere.  In  Figure  1  it  can  be  seen  that  at  low  temperature 
the  energy  of  the  cluster  is  Indeed  independent  of  the  size 
of  the  containment  sphere,  to  within  the  statistical  uncertainty 
of  the  calculation.  The  equilibrium  constants  are  shown  in 
Figure  3.  It  is  easy  to  show  from  Eqns.  (2.12)  and  (2.13) 
that,  at  each  p,  Kn(£)  is  an  increasing  function  of  the  radius 
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of  the  containment  sphere.  Thus,  in  Figure  3  the  dashed 
curve  should  always  be  above  the  solid  curve,  if  both  were 
calculated  exactly.  In  fact,  the  curves  cross,  but  the  dif¬ 
ference  between  the  two  curves  at  low  temperatures  is  approxi 
mately  equal  to  the  estimated  statistical  error  in  the  two 
calculations.  Thus,  to  within  the  statistical  uncertainty 
in  our  calculations,  the  equilibrium  constants  are  the  same 
for  both  choices  of  the  containment  sphere  radius. 


VI.  Results 


The  results  for  the  logarithms  of  the  equilibrium 
coefficients  are  displayed  in  Figure  4.  In  each  case,  these 
are  graphed  relative  to  In  Kn(0),  the  infinite  temperature 
result  [see  (2.17)]. 

The  equilibrium  coefficients  for  dimer  and  trimer  forma¬ 
tion  may  be  compared  indirectly  with  experiment  as  a  test  of 
the  method  and  of  the  intermolecular  potentials.  In  Section 
II  we  discuss  the  relationships  that  exist  between  the  second 
and  third  virial  coefficients,  I^Cp)  and  B^(p),  and  the 
second  and  third  equilibrium  coefficients,  K^Cp)  and  K^(P). 

These  allow  the  calculation  of  a  theoretical  from 

^(p),  as  well  as  the  calculation  of  an  experimental  In  K^(P) 
from  experimental  I^P)  and  B^(P). 

18*20 

Figure  5  provides  the  comparison  of  experimental 
and  theoretically  derived  second  virial  coefficients.  The 
experimental  values  are  shown  over  the  entire  range  of  tempera¬ 
tures  for  which  the  experiments  were  performed.  The  theo¬ 
retically  derived  curve  for  came  from  Eq.  (2.18)  plus  small 
corrections  for  the  effect  of  the  long-ranged  dipolar  part 
of  the  intermolecular  potential.  (The  second  virial  coeffi¬ 
cient  can  be  expressed  as  the  integral  of  the  Mayer  f  function 
over  all  relative  positions  and  orientations  of  two  molecules. 
Equation  (2.18)  for  B2  was  derived  assuming  that  the  interaction. 


and  hence  the  Mayer  f  function,  is  zero  for  intermolecular 
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distances  larger  than  2R2-  Thus,  to  obtain  B2  for  the  full 

Watts  potential,  it  is  necessary  to  add  the  contribution 

from  the  integral  for  distances  larger  than  2R2-  The  latter 

21 

was  calculated  analytically  assuming  that  the  Watts  poten¬ 
tial  was  dipolar  at  those  large  distances.) 

The  agreement  between  the  experimental  and  theoretical 
second  virial  coefficients  is  seen  to  be  quite  good.  Over 
the  temperature  range  for  which  experimental  information  is 
available,  the  differences  between  the  theory  and  experiment 
are  only  about  as  large  as  the  differences  among  the  various 
experiments.  It  should  be  noted,  however,  that  the  Watts 
potential  was  constructed  to  give  second  virial  coefficients 
in  agreement  with  experiment.  The  agreement  between  experi¬ 
ment  and  theory  as  shown  in  Figure  5  should  therefore  be 

interpreted  as  confirmation  of  our  method  for  computing  phase 

14 

integrals .  The  procedure  used  by  Watts  to  compute  second 
virial  coefficients  for  various  water  potentials  kept  the 
water  molecules  rigid  at  their  equilibrium  internuclear  con¬ 
figurations.  The  present  calculations  make  no  such  restric¬ 
tion,  and  so  the  present  method  for  calculating  dimer  forma¬ 
tion  constants  of  vibrating  molecules  gives  a  convenient  way 
of  calculating  second  virial  coefficients  of  vibrating  molecules. 

48 


1 


In  the  case  of  water,  however,  the  effect  of  vibrations  in 
this  purely  classical  calculation  of  the  second  virial  coef¬ 
ficient  is  small. 

Figure  6  compares  K^,  the  trimer  equilibrium  constant 
obtained  from  simulations  using  the  larger  than  standard 
trimer  containment  sphere,  with  an  experimental  value  of 
calculated  from  experimental  second  and  third  virial  coeffi¬ 
cients  using  Eq.  (2.20).  (In  order  to  be  consistent  with 
this  equation,  we  have  subtracted  from  the  experimental  second 
virial  coefficients  an  estimate  of  the  long  ranged  dipolar  con¬ 
tribution,  since  and  K'  were  calculated  for  a  potential 
that  had  no  such  long  ranged  part.  In  principle,  the  corres¬ 
ponding  correction  to  the  experimental  third  virial  coeffi¬ 
cient  should  also  be  made,  but  that  is  much  more  difficult 
and  is  not  necessary  since  the  experimental  third  virial 
coefficient  makes  only  a  small  contribution  to  the  experimental 
curve  in  this  figure.)  The  difference  between  the  experimental 
and  theoretical  results  is  of  about  the  same  size  as  the  dif¬ 
ference  between  different  experiments  and  the  statistical 
uncertainties  in  the  calculations. 

In  principle  the  third  virial  coefficient  can  be  calcu¬ 
lated  from  and  using  Eq.  (2.20).  Unfortunately,  this  is 
not  possible  in  practice  because  of  the  magnitude  of  the 


statistical  error  in  and  because  Eq.  (2.20)  gives  the  third 
virial  coefficient  as  the  small  difference  between  much  larger 
numbers . 


VII.  Discussion 


In  this  paper  we  have  presented  a  practical  procedure 
for  evaluating  the  equilibrium  constants  for  the  formation  of 
clusters  in  the  gas  phase  using  equilibrium  classical  statis¬ 
tical  mechanics  and  assumptions  about  the  intermolecular  inter¬ 
actions.  The  calculations  were  based  on  a  specific  model 
potential  for  water  and  on  the  assumption  that  the  energy  of  a 
collection  of  molecules  is  the  sum  of  the  interactions  between 
pairs  of  molecules .  The  Watts  intermolecular  potential  is 

known  to  be  significantly  different  from  the  true  water  poten- 
22 

tial.  The  assumption  of  pairwise  additivity  of  the  potential 
is  also  known  to  be  not  quantitatively  accurate  for  water. 
Moreover,  for  water  the  internal  motions  cannot  be  regarded  as 
classical.  Because  of  these  three  facts,  the  calculations 
presented  here  should  not  be  regarded  as  predictions  of  the 
equilibrium  constants  for  water;  rather  they  are  example  calcu¬ 
lations  designed  to  show  the  feasibility  of  the  method.  The 
method  could  easily  be  applied  to  a  different  choice  of  poten¬ 
tials,  pairwise  nonadditive  potentials,  clusters  of  ions  and 
molecules,  and  other  situations.  The  restriction  to  classical 
mechanics  is  intrinsic  to  the  molecular  dynamics  method. 

9 

The  method  is  an  alternative  to  the  methods  of  Lee  ejt  al. 
and  of  Mruzik  ejt  al.  ^  It  differs  from  that  of  Lee  ejt  al.  by 


51 


1 


using  an  integration  to  infinite  temperature  to  establish  the 
absolute  free  energy  of  a  cluster  rather  than  using  an  inte¬ 
gration  over  the  value  of  the  size  of  the  containment  sphere. 

It  differs  from  the  method  of  Mruzik  e_t  ajL.  in  that  the  latter 
calculates  the  free  energy  difference  for  clusters  differing 
by  one  molecule  and  relies  on  the  fact  that  the  absolute  free 
energy  of  the  single  molecule  cluster  can  be  evaluated  exactly. 

As  a  byproduct  of  this  method,  we  have  a  simulation 
method  for  calculating  the  second  virial  coefficient  of  non- 
rigid  molecules. 
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Appendix:  Velocity  Form  of  the  Verlet  Algorithm 


In  this  appendix  we  state  the  velocity  form  of  the  Verlet 
23 

algorithm,  which  we  use  for  integrating  Newton's  equations 
of  motion. 

Assuming,  for  notational  simplicity,  that  there  is  one 
degree  of  freedom,  Newton's  equations  can  be  written  as 

i:  =  f  (r) 

in  which  f(r)  is  the  force  divided  by  the  mass.  Let  h  be  the 
time  step  for  the  numerical  integration,  and  let 


t  *=  nh 
n 


n 


n 


*<*»> 


The  Verlet  algorithm  is  based  on  the  following  approximations 


in“  (rn+l‘rn-l)/2h 


£  **  (r  ,  i  -  2r  +  r  , 

n  v  n+1  n  n-1 

The  Verlet  algorithm 
as  an  equality  to  get 


)/h2. 

is  obtained  by  treating  Eq. 


(A.  1) 
(A.  2) 

(A.  2) 


n+1 


2r  -  r  ,  +h2f  (r  ). 
n  n-x  n 


(A. 3) 


This  can  be  the  basis  for  a  recursive  procedure  for  calcu¬ 
lating  all  subsequent  rR  from  Tq  and  r^.  If  the  velocity  is 


needed  we  define 


v 

n 


(r 


n+1 


V!)/2h 


(A.4) 


and  using  Eq .  (A.l)  we  interpret  as  r  .  In  the  numerical 

analysis  literature,  the  Verlet  algorithm  is  known  as  the 

24 

"explicit  central  difference  method." 

A  better  way  of  implementing  the  Verlet  algorithm  on  a 
computer  with  finite  precision  is  to  use  the  "summed  form." 
Defining 


z 

n 


rn>/h 


it  is  easy  to  show  that  Eq.  (A. 3)  is  equivalent  to 


r  +hz  i 
n-1  n-1 

z  .  +  hf  (r  )  . 
n-l  n 


These  equations  can  be  iterated  to  obtain  all  subsequent 

values  of  r  and  z  from  and  zn.  If  velocities  are  needed, 
n  n  0  0 

they  can  be  obtained  from 


v 

n 


(z  +z  ,)/2. 
n  n-i 


These  equations  are  mathematically  equivalent  to  the  Verlet 

algorithm,  but  they  are  not  numerically  equivalent  and  are 

o  / 

superior  on  a  computer  of  finite  precision. 
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The  velocity  form  of  the  Verlet  algorithm  is  one  in  which 


appears  directly  in  the  equations  to  be  iterated.  It  is 
straightforward  to  show  that  the  following  equations  are 
equivalent  to  Eqs .  (A. 3)  and  (A. A). 


r  ,  «  r  +  hv  +  n  f (r  )/2 
n+1  n  n  v  n 


v„+l  "  vn  +  h[f<rn+l)  +  f(rn)1/2- 


These  equations  retain  the  superior  numerical  precision  of 

the  summed  form.  Since  they  are  a  way  of  directly  getting  the 

position  and  velocity  at  the  end  of  the  time  step  from  the 

position  and  velocity  at  the  beginning  of  the  step,  they  provide 

an  easy  way  of  grafting  stochastic  collisions  into  the  algorithm. 

The  effect  of  a  stochastic  collision  is  merely  to  change  the 

value  of  v  just  before  r  is  to  be  calculated, 
n  n+1 
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Figure  Captions 


Figure  1 

Figure  2 


Figure  3 


Figure  U 
Figure  5 


U  -  nU,  for  n  *  2  to  n  *  5.  The  dashed  line  shows 
n  i 

U3  -  3U^  with  U3  calculated  using  the  larger  than 
standard  trimer  containment  sphere.  (Vertical  lines 
are  error  estimates  of  ±1  standard  deviation.) 

.  Dimer  (n  =  2)  and  pentamer  (n  =  5)  potential  energy 

heat  capacity,  -dU^/dp  (solid  lines),  and  fluctuation 

2  2 

in  potential  energy,  (U^  )  -  (U^)  (circles),  as  a 
function  of  p.  See  the  text  for  the  use  of  this 
graph  as  a  consistency  check  on  the  calculation. 

Ln[K3(p)/K3(0)]  calculated  using  both  the  standard 
containment  sphere  size  (solid  curve)  and  the  larger 
than  standard  containment  sphere  size  (dashed  curve). 
K3(0),  to  which  both  curves  are  referenced  is  that 
of  the  standard  trimer  containment  sphere. 

Ln[K  (p)/K  (0)]  for  n  *  2  to  n  *=  5. 
n  n 

The  second  virial  coefficient,  I^P),  from  experi¬ 
ment  and  theory.  The  solid  line  shows  calculated 
from  K^p).  The  dashed  line  is  from  the  data  of 
Keyes,  et  al.  (Reference  18);  the  dotted  line  is 
from  the  data  of  Kell  et  al.  (Reference  19);  the 


dot-dashed  line  is  from  the  data  of  Vukalovich 
et  al.  (Reference  20). 

Figure  6.  LnlK^fO/K^O)]  as  computed  from  theory  (solid 

curve)  and  from  experiments.  The  theoretical  curve 
was  obtained  from  molecular  dynamics  results  using 
the  larger  than  standard  trimer  containment  sphere. 
The  experimental  curves  were  obtained  using  Eq. 
(2.20)  and  second  and  third  virial  coefficient  data 
of  Kell  £t  a^l.  (Reference  19)  and  of  Vukalovich 
et  al .  (Reference  20) .  The  data  of  Kell  et  al . 
determine  the  dotted  curve  and  that  of  Vukalovich 
et  al.  determine  the  dot-dashed  curve. 
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